Project

rm(list=ls())
suppressPackageStartupMessages({
library(rjags)
library(bayesplot)
library(jagsUI)
library(glue)
library(loo)
library(knitr)
library(kableExtra)
  }
)

Bayesian Learning & Monte-Carlo Simulation final project

Introduction

Setting up environment, importing the dataset and visualising a graph of the two time series we will have to fit.

rm(list=ls())

# If set to TRUE higher order models will also be trained
intense = TRUE

# New Family Houses Sold: United States
# Source: https://fred.stlouisfed.org/series/HSN1F
data      <- read.csv("gdp_inflation.csv",header=T)
data$DATE <- as.Date(data$DATE)

# Plotting time series for visualisation purposes
plot(
  data$DATE[1:305],
  data$GDP_PC1[1:305],
  type="l",
  xlab="",
  ylab="",
  main="Percent change YoY (quarterly)"
)
lines(
  data$DATE[1:305],
  data$CPIAUCSL_PC1[1:305],
  type="l",
  col="red"
)

legend(x="bottom",legend=c("GDP","CPI"),col=c("black","red"),lty=1)

At this point, we setup our environment by converting the values in both time series into numericals, and then splitting them into training and testing. Finally, we define some useful variables which we will be using through out our script.

# Dataset extraction
x_axis      <- data$DATE[1:305]
series_GDP  <- as.numeric(data$GDP_PC1[1:305])
series_CPI  <- as.numeric(data$CPIAUCSL_PC1[1:305])

# Parameter setup (RMK: we assume both time series of same length)
TRAIN_PERC  <- 0.15
N_tot       <- length(series_GDP)
N_test      <- floor(TRAIN_PERC * N_tot)
N_train     <- N_tot - N_test

To further inspect our data we can also plot the respective auto-correlation functions:

acf(series_GDP, lag.max=100)

acf(series_CPI, lag.max=100)

The auto-correlation functions of the two proceses do not clearly match any given noticeable pattern, hence a more sophisticated model selection will be needed, where multiple models will be tested and evaluated one against the other.

Utility functions

For our convenience, we define a few utility functions which will be useful later. The first one automatises the interactions with JAGS, the second one will draw plots of both the whole time-series and its out-of-sample portion only, the third one will specifically plot a time-series along with the anomalies which have been detected through out its history, and so on.

fit_JAGS <- function(
  model_string,   # JAGS model string
  data_list,      # List of input values for JAGS
  to_save,        # Vector with parameter names which should be saved
  n_iter=5000,    # Number of iterations. Should be greater than n_burnin
  n_adapt=1000,   # JAGS n.adapt internal parameter
  n_chain=1,      # Number of chains to run
  n_burnin=1000,  # Burn-in value
  n_thin=1,       # Thinning value
  verbose=FALSE   # Displays JAGS logging / output (turn on for debugging)
) {
  
  output <- jags(
    model.file=textConnection(model_string),
    data=data_list,
    parameters.to.save=to_save,
    n.adapt=n_adapt,
    n.iter=n_iter,
    n.chains=n_chain,
    n.burnin=n_burnin,
    n.thin=n_thin,
    verbose=verbose
  )
  
  return(output)
}

visualise <- function(
  x_axis,           # Values to be displayed on the x-axis of the plots
  original_series,  # Original time series with all samples
  output,           # Output from JAGS, with in-sample and out-of-sample predictions
  title="",         # Title to display on the graphs
  keep=5            # Training sampled to keep in the test graph
) {
  # Length of total series
  N <- length(x_axis)
  
  
  # Plotting entire graph (training-set + predictions)
  # In-sample predicted values with their 2.5%-97.5% quantiles
  yp_in_mean  <- output$mean$yp_in
  yp_in_q1    <- output$q2.5$yp_in
  yp_in_q2    <- output$q97.5$yp_in
  in_length   <- length(yp_in_mean)
  
  # Out-of-sample predicted values with their 2.5%-97.5% quantiles and 25%-75% quantiles
  yp_out_mean <- output$mean$yp_out
  yp_out_q1   <- output$q2.5$yp_out
  yp_out_q2   <- output$q97.5$yp_out
  yp_out_q3   <- output$q25$yp_out
  yp_out_q4   <- output$q75$yp_out
  
  # Removing potentially non-finite values from quantiles
  yp_in_q1[is.infinite(yp_in_q1)]   <- NA
  yp_in_q2[is.infinite(yp_in_q2)]   <- NA
  yp_out_q1[is.infinite(yp_out_q1)] <- NA
  yp_out_q2[is.infinite(yp_out_q2)] <- NA
  yp_out_q3[is.infinite(yp_out_q3)] <- NA
  yp_out_q4[is.infinite(yp_out_q4)] <- NA
  
  # Computing extreme values for first graph
  min_y <- min(original_series, yp_in_q1, yp_out_q1)
  max_y <- max(original_series, yp_in_q2, yp_out_q2)
  
  # Plotting REAL time series (as points only), both training set and test-set
  plot(
    x_axis,
    original_series,
    main=title, col="red",
    xlab="time", ylab="price",
    ylim=c(min_y, max_y)
  )
  lines(
    x_axis,
    original_series,
    col="red"
  )
  
  # Plotting a vertical line signalling end of training set and horizontal
  # Showing the computer mean
  abline(v=x_axis[in_length], col="magenta", lwd=1, lty=2)
  abline(h=output$mean$m0,    col="green",   lwd=1, lty=2)
  
  # Plotting in-sample predictions as well as 2.5%-97.5% quantiles
  points(x_axis[1:in_length], yp_in_mean, col="blue", pch="*")
  lines(x_axis[1:in_length],  yp_in_q1,   col="blue", lwd=1.5)
  lines(x_axis[1:in_length],  yp_in_q2,   col="blue", lwd=1.5)
  
  # Plotting out-of-sample predicted values as well as 2.5%-97.5%
  # quantiles and 25%-75% quantiles
  points(x_axis[(in_length+1):N], yp_out_mean, col="orange", pch="*")
  lines(x_axis[(in_length+1):N],  yp_out_q1,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q2,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q3,   col="orange", lwd=1.5, lty=2)
  lines(x_axis[(in_length+1):N],  yp_out_q4,   col="orange", lwd=1.5, lty=2)
  

  # Plotting onyl out-of-sample portion of graph (with last 5 training samples)
  # Re-computing extreme values for second graph
  min_y <- min(original_series[(in_length+1 - keep):N], yp_out_q1)
  max_y <- max(original_series[(in_length+1 - keep):N], yp_out_q2)
  
  # Plotting REAL time series (as points only), both training set and test-set
  plot(
    x_axis[(in_length+1 - keep):N],
    original_series[(in_length+1 - keep):N],
    main=title, col="red",
    xlab="time", ylab="price",
    ylim=c(min_y, max_y)
  )
  lines(
    x_axis[(in_length+1-keep):N],
    original_series[(in_length+1-keep):N],
    col="red"
  )
  
  # Plotting a vertical line signalling end of training set
  abline(v=x_axis[in_length+1], col="magenta", lwd=1, lty=2)
  abline(h=output$mean$m0,      col="green",   lwd=1, lty=2)
  
  # Plotting out-of-sample predicted values as well as 2.5%-97.5%
  # quantiles and 25%-75% quantiles
  points(x_axis[(in_length+1):N], yp_out_mean, col="orange", pch="*")
  lines(x_axis[(in_length+1):N],  yp_out_q1,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q2,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q3,   col="orange", lwd=1.5, lty=2)
  lines(x_axis[(in_length+1):N],  yp_out_q4,   col="orange", lwd=1.5, lty=2)
}

# Function in charge of visualising VAR models (requires custom rendering).
# This acts as a wrapper function to the standard visualise(...) function, defined above.
visualise_VAR = function(
  x_axis,               # Values to be displayed on the x-axis of the plots
  original_series_gdp,  # Original GDP time series with all samples
  original_series_infl, # Original INFLATION time series with all samples
  output,               # Output from VAR JAGS, with in-sample and out-of-sample predictions
  title="VAR",          # Title to display on the graphs
  keep=5                # Training sampled to keep in the test graph
) {
  
  # Marginalise GDP
  output_gdp = list()
  output_gdp$mean$yp_in   = output$mean$yp_in_gdp
  output_gdp$q2.5$yp_in   = output$q2.5$yp_in_gdp
  output_gdp$q97.5$yp_in  = output$q97.5$yp_in_gdp
  output_gdp$mean$yp_out  = output$mean$yp_out_gdp
  output_gdp$q2.5$yp_out  = output$q2.5$yp_out_gdp
  output_gdp$q97.5$yp_out = output$q97.5$yp_out_gdp
  output_gdp$q25$yp_out   = output$q25$yp_out_gdp
  output_gdp$q75$yp_out   = output$q75$yp_out_gdp
  output_gdp$mean$m0      = output$mean$m_gdp
  
  # Marginalise CPI
  output_infl = list()
  output_infl$mean$yp_in   = output$mean$yp_in_infl
  output_infl$q2.5$yp_in   = output$q2.5$yp_in_infl
  output_infl$q97.5$yp_in  = output$q97.5$yp_in_infl
  output_infl$mean$yp_out  = output$mean$yp_out_infl
  output_infl$q2.5$yp_out  = output$q2.5$yp_out_infl
  output_infl$q97.5$yp_out = output$q97.5$yp_out_infl
  output_infl$q25$yp_out   = output$q25$yp_out_infl
  output_infl$q75$yp_out   = output$q75$yp_out_infl
  output_infl$mean$m0      = output$mean$m_infl
  
  # Call visualisation on marginalisations
  visualise(
    x_axis,
    original_series_gdp,
    output_gdp,
    title=glue("GDP series, {title}"),
    keep
  )
  visualise(
    x_axis,
    original_series_infl,
    output_infl,
    title=glue("CPI series, {title}"),
    keep
  )
}

# Function in charge of visualising the anomaly points on top of its corresponding time series
visualise_anomalies <- function(
  x_axis,             # X-axis to use for the plot
  original_series,    # Original time series to display
  high_var_indices,   # List of indices corresponding to high variance points
  med_var_indices,    # List of indices corresponding to medium variance points
  title="",           # Name of the time series used
  start_point=1       # Value after which plotting should start (in order to obtain close up graphs)
) {
  # Plot actual time series
  N <- length(x_axis)
  plot(
    x_axis[start:N],
    original_series[start:N],
    main=glue("Anomalies for {title}"),
    t="l", xlab="time", ylab="price"
  )

  # Plotting strong anomaly points
  points(x_axis[high_var_indices], original_series[high_var_indices], col="red")
  for (high_anomaly in high_var_indices) {
    abline(v=x_axis[high_anomaly], col="red", lwd=1.5)
  }
  
  # Plotting medium anomaly points
  points(x_axis[med_var_indices], original_series[med_var_indices], col="orange")
  for (medium_anomaly in med_var_indices) {
    abline(v=x_axis[medium_anomaly], col="orange", lty=3, lwd=1.5)
  }
}

# Function in charge of classifying anomaly points given their delta value after running the anomaly detection simulation
classify_anomalies <- function(
  output,             # JAGS simulation output
  N_tot,              # Total number of samples
  title="",           # Name of time series to display
  high_value_thr=0.2, # High variance threshold value used
  med_value_thr=0.60  # Medium variance threshold value used
) {
  # Obtaining precision, variance and delta
  tau     <- output$mean$tau
  sigma2  <- 1 / tau
  delta   <- output$mean$d
  
  # Classifying the values of d
  high_var_indices  <- seq(1:(N_tot-1))[delta < high_value_thr]
  med_var_indices   <- seq(1:(N_tot-1))[high_value_thr < delta & delta < med_value_thr]
  low_var_indices   <- seq(1:(N_tot-1))[delta > med_value_thr]
  
  # Logging counts
  n_high  <- length(high_var_indices)
  n_med   <- length(med_var_indices)
  print(
    glue("Found {n_high} strong and {n_med} medium outliers out {N_tot} total samples for {title}")
  )
  
  # Returning the indices
  return(list(
    high_var_indices = high_var_indices,
    med_var_indices  = med_var_indices,
    low_var_indices  = low_var_indices
  ))
}

# Function in charge of plotting the anomaly point delta values with their relative threshold lines
display_classification <- function(
  x_axis,             # X-axis to use for the plot
  delta,              # Sequence of delta values
  high_var_indices,   # List of indices corresponding to high variance points
  med_var_indices,    # List of indices corresponding to medium variance points
  low_var_indices,    # List of indices corresponding to low variance points
  title="",           # Name of the time series used
  high_value_thr=0.2, # High variance threshold value used
  med_value_thr=0.60  # Medium variance threshold value used
) {
  # Plotting the posterior values of delta and classifying them
  plot(
    x_axis[low_var_indices], delta[low_var_indices],
    xlab="time",
    ylab="Posterior delta values",
    main=glue("{title} anomaly samples"),
    xlim=c(min(x_axis), max(x_axis)),
    ylim=c(0, 1)
  )
  
  points(x_axis[high_var_indices], delta[high_var_indices], col="red")
  points(x_axis[med_var_indices], delta[med_var_indices]  , col="orange")
  abline(h=high_value_thr, col="red")
  abline(h=med_value_thr, col="orange")
}

# K is the number of effective params
compute_metrics <- function(output, series, N_train, k) {
  
  ############################################################
  ###################### IN-SAMPLE ###########################
  ############################################################
  # In-sample: DIC, WAIC, BIC
  ll_mat <- output$sims.list$loglik
  DIC_val <- compute_dic_marginal(ll_mat)
  
  waic_res <- suppressWarnings( waic(ll_mat) )
  WAIC_val <- waic_res$estimates["waic","Estimate"]
  
  #BIC computation
  #Sum of the log-likelihoods for each draw
  ll_sum <- rowSums(ll_mat)
  
  # maximum of the sum: estimate of log p(y|theta_hat), that is, log-lik_max
  lhat <- max(ll_sum)
  
  # BIC_val <- -2 * lhat + k * log(N_train)
  BIC_val <- -2 * lhat + k * log(N_train - output$mean$train_metrics_skip)
  
  # BIC Approx
  # BIC_val  <- DIC_val + k * (log(N_train) - 2) 
  
  ############################################################
  #################### OUT-OF-SAMPLE #########################
  ############################################################

  y_pred <- if (is.null(output$mean$yp_onestep_out)) output$mean$yp_out else output$mean$yp_onestep_out

  y_true <- series[(N_train+1):length(series)]
  e <- y_pred - y_true
  
  # MSE
  MSE_val <- mean((e)^2)
  
  #MAE
  MAE_val <- mean(abs(e))
  MAPE_val <- mean(abs(e / y_true)) * 100
  
  #SMAPE
  SMAPE_val <- mean(2 * abs(e) / (abs(y_true) + abs(y_pred))) * 100
  
  #MASE
  scale <- mean(abs(diff(series[1:N_train])))
  MASE_val <- MAE_val / scale
  
 list(
    DIC   = round(DIC_val, 1),
    WAIC  = round(WAIC_val, 1),
    BIC   = round(BIC_val, 1),
    MSE   = round(MSE_val, 3),
    MAE   = round(MAE_val, 3),
    MAPE  = round(MAPE_val, 1),
    SMAPE = round(SMAPE_val, 1),
    MASE  = round(MASE_val, 3)
  )
}

compute_dic_marginal <- function(ll_mat) {
  dev      <- -2 * ll_mat
  dev_sum  <- rowSums(dev)
  Dbar     <- mean(dev_sum)
  pD       <- 0.5 * var(dev_sum)
  DIC      <- Dbar + pD
  return(DIC)
}

# Appends metrics row to provided metrics_df and returns it
# If no arguments are passed, generates and returns a new empty metrics_df
metrics_df_append <- function(metrics_df = NA, new_mets = NA, name="...") {
  if (anyNA(metrics_df)) {
    return(data.frame(
      Model = character(),
      DIC   = numeric(),  
      WAIC  = numeric(), 
      BIC   = numeric(),
      MSE   = numeric(),  
      MAE   = numeric(), 
      MAPE  = numeric(),
      SMAPE = numeric(),  
      MASE  = numeric(),
      stringsAsFactors = FALSE
    ))
    
    
  } else if (!anyNA(metrics_df) && !all(is.na(new_mets))) {
    
    return(rbind(
      metrics_df,
      data.frame(
        Model = name,
        DIC   = new_mets$DIC,
        WAIC  = new_mets$WAIC,
        BIC   = new_mets$BIC,
        MSE   = new_mets$MSE,
        MAE   = new_mets$MAE,
        MAPE  = new_mets$MAPE,
        SMAPE = new_mets$SMAPE,
        MASE  = new_mets$MASE,
        stringsAsFactors = FALSE
      ))
    )
  }
  
  return(metrics_df)
}

JAGS model strings

At this point, we define all our JAGS model strings so as to be able to reuse them multiple times through out the script. We start with with defining the model for a generic \(\text{AR}(n)\) sequence. The auto-regressive order \(n\) is passed to the model by means of the order variable.

\[ y(t) = \left[ \sum_{i=1}^{n} \alpha_i y(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2) \]

For any order \(n ≥ 1\) the following model will correctly implement an \(\text{AR}(n)\) process, which allows us to automatically test multiple AR processes easily and at once.

# AR(n) string
modelAR_string <- "
model {
  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################
  
  # 1. Bootstrapping the AR process
  mu[1:order]     <- y[1:order]
  yp_in[1:order]  <- y[1:order]
  
  # 2. Recursive step of the AR process
  for (t in (order+1):N_train) {
    mu[t]     <- inprod(alpha_rev, y[(t-order):(t-1)]) + m0
    y[t]      ~  dnorm(mu[t], tau)
    yp_in[t]  ~  dnorm(mu[t], tau)
  }
  
  for (t in (train_metrics_skip+1):N_train) {
    loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
  }
  
  ############################################################
  ######################## PREDICTION ########################
  ############################################################
  
  # Init helper array z
  z[1:order] <- y[(N_train-order+1):N_train]
  
  for (t in (order+1):(order+N_test)) {
    mu_z[t] <- inprod(alpha_rev, z[(t-order):(t-1)]) + m0
    z[t]    ~  dnorm(mu_z[t], tau)
  }
  
  yp_out <- z[(order+1):(order+N_test)]
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  for (t in 1:N_test) {
    mu_onestep[t]     <- inprod(alpha_rev, y[(N_train+t-order):(N_train+t-1)]) + m0
    yp_onestep_out[t]  ~ dnorm(mu_onestep[t], tau)
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  m0  ~ dnorm(0.0, 1.0E-4)
  tau ~ dgamma(0.1, 10)
  for (i in 1:order) {
    alpha[i] ~ dunif(-1.5, 1.5)
    alpha_rev[order - i + 1] <- alpha[i]
  }
  
  # Defining converstion from precision to variance here
  sigma2 <- 1 / tau
}"

We now define a generic \(\text{MA}(n)\) process. The memory order \(n\) is passed to the model by means of the order variable.

\[ y(t) = \left[ \sum_{i=1}^n \alpha_i \eta(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2) \]

For any order \(n ≥ 1\) the following model will correctly implement a \(\text{MA}(n)\) process, which allows us to automatically test multiple MA processes easily and at once.

modelMA_string <- "
model {

  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################

  e[1:order] = rep(0, order)
  yp_in[1:order]  <- y[1:order]
  
  # 2. Recursive step of the MA process
  for (t in (order+1):N_train) {
    mu[t]    <- m0 + inprod(beta_rev, e[(t-order):(t-1)])
    y[t]      ~ dnorm(mu[t], tau)
    yp_in[t]  ~ dnorm(mu[t], tau)
    e[t]     <- y[t] - mu[t]
  }

  for (t in (train_metrics_skip+1):N_train) {
    loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
  }
  
  ############################################################
  ######################## PREDICTION ########################
  ############################################################

  for (t in 1:N_test) {
    e[N_train+t] ~ dnorm(0, tau)
    mu_out[t]   <- m0 + inprod(beta_rev, e[(N_train+t-order):(N_train+t-1)])
    yp_out[t]    ~ dnorm(mu_out[t], tau)
  }
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  e_onestep[1:order]   <- e[(N_train-order+1):N_train]
  for (t in 1:N_test) {
    mu_onestep[t]      <- m0 + inprod(beta_rev, e_onestep[t:(t+order-1)])
    yp_onestep_out[t]   ~ dnorm(mu_onestep[t], tau)
    e_onestep[t+order] <- y[N_train+t] - mu_onestep[t]
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  m0    ~ dnorm(0.0, 1.0E-4)
  tau   ~ dgamma(0.01, 0.01)

  for (i in 1:order){
    beta[i]      ~ dnorm(0, 1/10)
    beta_rev[i] <- beta[order-i+1]
  }
  
  # Defining converstion from precision to variance here
  sigma2 <- 1 / tau
}"

Finally, the following model will merge the two previous models defining an \(\text{ARMA}(p, q)\) model.

\[ y(t) = \left[ \sum_{i=1}^p \beta_i y(t-i) \right] + \left[ \sum_{i=1}^q \alpha_i \eta(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2). \]

The parameters order_ar and order_ma allow us to specify the two orders when simulating the model so as to be able to test all possible combinations of auto-regressive and moving-average processes.

modelARMA_string <- "
model {
  
  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################

  e[1:order]     <- rep(0, order)
  mu[1:order]    <- y[1:order]
  yp_in[1:order] <- y[1:order]
  
  # 2. Recursive step of the MA process
  for (t in (order+1):N_train) {
    mu_ar[t] <- inprod(alpha_rev, y[(t-order):(t-1)])
    mu_ma[t] <- inprod(beta_rev, e[(t-order):(t-1)])
    mu[t]    <- m0 + mu_ar[t] + mu_ma[t]
    y[t]      ~ dnorm(mu[t], tau)
    yp_in[t]  ~ dnorm(mu[t], tau)
    e[t]     <- y[t] - mu[t]
  }

  for (t in (train_metrics_skip+1):N_train) {
    loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
  }
  
  ############################################################
  ######################## PREDICTION ########################
  ############################################################

  z[1:order]   <- y[(N_train-order+1):N_train]
  e_z[1:order] <- e[(N_train-order+1):N_train]

  for (t in (order+1):(order+N_test)) {
    e_z[t]        ~ dnorm(0, tau)
    mu_ar_out[t] <- inprod(alpha_rev, z[(t-order):(t-1)])
    mu_ma_out[t] <- inprod(beta_rev, e_z[(t-order):(t-1)])
    mu_z[t]      <- m0 + mu_ar_out[t] + mu_ma_out[t]
    z[t]          ~ dnorm(mu_z[t], tau)
  }

  yp_out <- z[(order+1):(order+N_test)]
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  e_onestep[1:order]   <- e[(N_train-order+1):N_train]
  for (t in 1:N_test) {
    mu_ar_onestep[t] <- inprod(alpha_rev, y[(N_train+t-order):(N_train+t-1)])
    mu_ma_onestep[t]   <- inprod(beta_rev, e_onestep[t:(t+order-1)])
    mu_onestep[t]      <- m0 + mu_ar_onestep[t] + mu_ma_onestep[t]
    yp_onestep_out[t]   ~ dnorm(mu_onestep[t], tau)
    e_onestep[t+order] <- y[N_train+t] - mu_onestep[t]
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  m0    ~ dnorm(0.0, 1.0E-4)
  tau   ~ dgamma(0.01, 0.01)
  
  for (i in 1:order_ar) {
    alpha[i] ~ dunif(-1.5, 1.5)
  }
  for (i in (order_ar+1):order) {
    alpha[i]      <- 0
  }
  for (i in 1:order) {
    alpha_rev[i] <- alpha[order-i+1]
  }

  for (i in 1:order_ma) {
    beta[i] ~ dnorm(0, 1/10)
  }
  for (i in (order_ma+1):order) {
    beta[i]      <- 0
  }
  for (i in 1:order) {
    beta_rev[i] <- beta[order-i+1]
  }
  
  # Defining converstion from precision to variance here
  sigma2 <- 1 / tau
  
}"

Modelling

At this point we can start modelling our two time series using the stochastic processes for which we defined the previous model strings. We will start with the GDP time-series and then focus on the CPI time-series. Finally, we will also attempt to model a correlated process using a \(\text{VAR}(n)\) model, where both time-series will be modelled together at once.

GDP time series

# Metrics dataframe setup
metrics_df_GDP <- metrics_df_append()

We will, at first, start off with some \(\text{AR}(n)\) models, with increasing values of \(n\). We will use the obtained results ot decide whether or not the selected model may apply to the problem at hand, and we will also evaluate which regressive order is best.

to_save <- c(
  "alpha",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order <- if (intense) 9 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_GDP,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelAR_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_GDP, output, title=glue("GDP series, AR({order})"))
  
  # Update metrics
  mets <- compute_metrics(output, series_GDP, N_train, k=order+2)
  metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("AR({order})"))
}

We will now also fit some \(\text{MA}(n)\) processes (for now \(1 \leq n \leq 5\)) to evaluate whether anyone might be better.

to_save <- c(
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our MA(n) from n=1 to max_order
max_order <- if (intense) 5 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_GDP,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelMA_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_GDP, output, title=glue("GDP series, MA({order})"))
  
  # Updating metrics
  mets <- compute_metrics(output, series_GDP, N_train, k=order+2)
  metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("MA({order})"))
}

Finally, we combine the auto-regressive and moving-average proccesses and try to model our series using an \(\text{ARMA}(p, q)\) process. For now, we will attempt modelling for all combinations of \(1 \leq p \leq 4\) and \(1 \leq q \leq 4\).

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order_ar <- if (intense) 4 else 2
max_order_ma <- if (intense) 4 else 2
for (order_ma in 1:max_order_ma) {
  for (order_ar in 1:max_order_ar) {
    data_list <- list(
      "y"=series_GDP,
      "N_train"=N_train,
      "N_test"=N_test,
      "order"=max(order_ar, order_ma),
      "order_ar"=order_ar,
      "order_ma"=order_ma,
      "train_metrics_skip"=50
    )
    
    # Perform simulations and visualise output
    output <- fit_JAGS(
      modelARMA_string,
      data_list,
      to_save
    )
    visualise(x_axis, series_GDP, output, title=glue("GDP series, ARMA({order_ar}, {order_ma})"))
    
    # Updating metrics
    mets <- compute_metrics(output, series_GDP, N_train, k=order_ar+order_ma+2)
    metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("ARMA({order_ar}, {order_ma})"))
  }
}

CPI time series

# Metrics dataframe setup
metrics_df_CPI <- metrics_df_append()

We continue our modelling moving to the CPI time-series. We will proceede as for the GDP, trying our a few \(\text{AR}(n)\) models, then moving to some \(\text{AR}(n)\) and finally, putting both together, some \(\text{ARMA}(p, q)\) models.

to_save <- c(
  "alpha",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order <- if (intense) 9 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_CPI,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelAR_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_CPI, output, title=glue("CPI series, AR({order})"))
  
  if (order == 9) {
    output_ar9_cpi = output
  }
  
  # Update metrics
  mets <- compute_metrics(output, series_CPI, N_train, k=order+2)
  metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("AR(",order,")"))
}

We will now also fit some \({MA}(n)\) processes (for now \(1 \leq n \leq 5\)) to evaluate whether anyone might be better.

to_save <- c(
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our MA(n) from n=1 to max_order
max_order <- if (intense) 5 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_CPI,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelMA_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_CPI, output, title=glue("CPI series, MA({order})"))
  
  # Updating metrics
  mets <- compute_metrics(output, series_CPI, N_train, k=order+2)
  metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("MA(",order,")"))
}

Finally, we combine the auto-regressive and moving-average proccesses and try to model our series using an \(\text{ARMA}(p, q)\) process. For now, we will attempt modelling for all combinations of \(1 \leq p \leq 4\) and \(1 \leq q \leq 4\).

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order_ar <- if (intense) 4 else 2
max_order_ma <- if (intense) 4 else 2
for (order_ma in 1:max_order_ma) {
  for (order_ar in 1:max_order_ar) {
    data_list <- list(
      "y"=series_CPI,
      "N_train"=N_train,
      "N_test"=N_test,
      "order"=max(order_ar, order_ma),
      "order_ar"=order_ar,
      "order_ma"=order_ma,
      "train_metrics_skip"=50
    )
    
    # Perform simulations and visualise output
    output <- fit_JAGS(
      modelARMA_string,
      data_list,
      to_save
    )
    visualise(x_axis, series_CPI, output, title=glue("CPI series, ARMA({order_ar}, {order_ma})"))
    
    # Updating metrics
    mets <- compute_metrics(output, series_CPI, N_train, k=order_ar+order_ma+2)
    metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("ARMA({order_ar}, {order_ma})"))
  }
}

VAR modelling

Finally, we will try to fit generic order \(\text{VAR}(n)\) models, which have the following form:

\[ \vec{y}(t) = \vec{\mu}_0 + \left[ \sum_{i=1}^n A_i \vec{y}(t-i) \right] + \vec{\eta}(t) \qquad \forall i \; A_i \in \mathbb{M}(2 \times 2, \mathbb{R}) \qquad \vec{\eta} \sim \mathcal{N}(\vec{0}, \Sigma) \]

Where \(\Sigma\) is the covariance matrix of our white noise (now an \(\mathbb{R}^2\) vector) and \(\vec{y}(t) \in \mathbb{R}^2\) now contains both the GDP and CPI time-series, which will be modelled together in a correlated fashion. By passing the parameter order during simulation we can specify the auto-regressive order \(n\) of the multivariate process.

The idea behind this model is to attempt to model not only the two time-series independently, but also in a correlated way, as we work under the assumption that the two may influence each other.

# Prior fragment for Wishart distributions (positive definite covariance)
prior_wishart <- "
  # Wishart prior for a full covariance matrix
  R[1,1] <- 1
  R[2,2] <- 1
  R[1,2] <- 0
  R[2,1] <- 0
  tau ~ dwish(R, 3) # df=3 is the minimum for a 2x2 matrix
"

# Prior fragment for independent Gamma distributions (diagonal covariance)
prior_diagonal <- "
  # Independent Gamma priors for diagonal elements of precision matrix
  tau[1,1] ~ dgamma(0.01, 0.01)
  tau[2,2] ~ dgamma(0.01, 0.01)
  tau[1,2] <- 0 
  tau[2,1] <- 0
"

modelVAR_string = "
model {

  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################
  
  mu[1:2,1:order] = y[,1:order]
  yp_in[1:2,1:order] = y[,1:order]
  
  for (t in (order+1):N_train) {
    mu_tmp[1:2, t, 1] <- m0
    for (j in 1:order) {
      mu_tmp[1:2,t,j+1] <- mu_tmp[1:2, t, j] + A[,,j] %*% y[,(t-j)]
    }
    mu[1:2,t]   <- mu_tmp[1:2,t,order+1]
    y[1:2,t]     ~ dmnorm(mu[1:2,t], tau)
    yp_in[1:2,t] ~ dmnorm(mu[1:2,t], tau)
  }

  for (t in (train_metrics_skip+1):N_train) {
    loglik_gdp[t-train_metrics_skip]  <- logdensity.norm(y[1,t], mu[1,t], 1 / sigma2[1,1])
    loglik_infl[t-train_metrics_skip] <- logdensity.norm(y[2,t], mu[2,t], 1 / sigma2[2,2])
  }

  ############################################################
  ######################## PREDICTION ########################
  ############################################################
  
  z[1:2,1:order] = y[,(N_train-order+1):N_train]
  
  for (t in (order+1):(order+N_test)) {
    mu_z[1:2, t, 1] <- m0
    for (j in 1:order) {
      mu_z[1:2,t,j+1] <- mu_z[1:2, t, j] + A[,,j] %*% z[,(t-j)]
    }
    z[1:2,t] ~ dmnorm(mu_z[1:2,t,order+1], tau)
  }
  
  yp_out[1:2,1:N_test] = z[1:2,(order+1):(order+N_test)]
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  for (t in 1:N_test) {
    mu_onestep_tmp[1:2, t, 1] <- m0
    for (j in 1:order) {
      mu_onestep_tmp[1:2,t,j+1] <- mu_onestep_tmp[1:2, t, j] + A[,,j] %*% y[,(N_train+t-j)]
    }
    mu_onestep[1:2,t]     <- mu_onestep_tmp[1:2,t,order+1]
    yp_onestep_out[1:2,t]  ~ dmnorm(mu_onestep[1:2,t], tau)
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################

  for (i in 1:2) {
    for (j in 1:2) {
      for (o in 1:order) {
        A[i, j, o] ~ dunif(-1.2, 1.2)
        # A[i, j, o] ~ dnorm(0, 10)
      }
    }
  }
  
  # Dynamically inserted prior for tau
  
  {{prior_for_tau}}
  
  # __if UNCORRELATED_ERRORS__
  # 
  # tau[1,1] ~ dgamma(0.01, 0.01)
  # tau[2,2] ~ dgamma(0.01, 0.01)
  # tau[1,2] <- 0
  # tau[2,1] <- 0
  
  # __else__
  # 
  # tau     ~ dwish(R, 3) # dwish(R=ScaleMatrix... Identity, df=degrees of freedom)
  
  m_gdp   ~ dnorm(0.0, 1.0E-4)
  m_infl  ~ dnorm(0.0, 1.0E-4)
  m0     <- c(m_gdp, m_infl)
  
  sigma2 <- inverse(tau)
  
  ############################################################
  ######################## EXTRACTION ########################
  ############################################################

  # Expose them for output in R
  loglik_gdp_out   <- loglik_gdp
  loglik_infl_out  <- loglik_infl

  yp_in_gdp   <- yp_in[1,]
  yp_in_infl  <- yp_in[2,]
  yp_out_gdp  <- yp_out[1,]
  yp_out_infl <- yp_out[2,]
  yp_onestep_out_gdp  <- yp_onestep_out[1,]
  yp_onestep_out_infl <- yp_onestep_out[2,]
}"

# Setting up parameters
y_var <- matrix(c(series_GDP, series_CPI), nrow=2, byrow=TRUE)
order <- 1
use_wishart <- TRUE

if (use_wishart) {
  prior_for_tau = prior_wishart
} else {
  prior_for_tau = prior_diagonal
}
final_modelVAR_string <- glue(modelVAR_string, .open="{{", .close="}}")
# Uncomment for debugging inspection
#print(glue(modelVAR_string, .open="{{", .close="}}"))

# Defining JAGS parameters
to_save <- c(
  "A",
  "m_gdp",
  "m_infl",
  "sigma2",
  "yp_in_gdp",
  "yp_in_infl",
  "yp_out_gdp",
  "yp_out_infl",
  "yp_onestep_out_gdp",
  "yp_onestep_out_infl",
  "train_metrics_skip",
  "loglik_gdp",   # marginal log-likelihood for GDP
  "loglik_infl"   # marginal log-likelihood for CPI
)
data_list <- list(
  "y"=y_var,
  "order"=order,
  "N_train" = N_train,
  "N_test" = N_test,
  "train_metrics_skip"=50
)

# Fitting VAR model and visualising it
jags_output_VAR = fit_JAGS(
  final_modelVAR_string,
  data_list,
  to_save,
  n_iter=6000,
  n_adapt=1000,
  n_chain=1,
  n_burnin=2000
)
visualise_VAR(
  x_axis,
  series_GDP,
  series_CPI,
  jags_output_VAR,
  title=glue("VAR({order})"),
  keep=5
)

# Wrapper for GDP-only metrics
output_gdp <- list(
  mean      = list(
    yp_out = jags_output_VAR$mean$yp_out_gdp,
    yp_onestep_out = jags_output_VAR$mean$yp_onestep_out_gdp,
    train_metrics_skip=jags_output_VAR$mean$train_metrics_skip
  ),
  sims.list = list(loglik  = jags_output_VAR$sims.list$loglik_gdp)
)

# Wrapper for CPI-only metrics
output_infl <- list(
  mean      = list(
    yp_out = jags_output_VAR$mean$yp_out_infl,
    yp_onestep_out = jags_output_VAR$mean$yp_onestep_out_infl,
    train_metrics_skip=jags_output_VAR$mean$train_metrics_skip
  ),
  sims.list = list(loglik  = jags_output_VAR$sims.list$loglik_infl)
)

# Compute univariate metrics for GDP and CPI as well
# RMK: just for visualisation purposes: can't compare marginalised metrics with univariate metrics
mets_VAR_gdp  <- compute_metrics(output_gdp,  series_GDP, N_train, k = 4*order + 2 + 3)
mets_VAR_infl <- compute_metrics(output_infl, series_CPI, N_train, k = 4*order + 2 + 3)

metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets_VAR_gdp, glue("VAR_univar_GDP({order})"))
metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets_VAR_infl, glue("VAR_univar_GDP({order})"))

VAR: Correlated vs Uncorrelated errors

We will now try to fit the \(\text{VAR}(1)\) model using two different priors for the precision matrix of \(\tau\) of the 2 errors.

Prior 1: diagonal precision matrix

We will force the errors to be uncorrelated. We will do this by setting the priors for the off-diagonal elements to be 0, while leaving a gamma prior to the diagonal elements, akin to the priors used of the univariate AR models.

Prior 2: positive definite precision matrix

We will set a Wishart prior to the precision matrix of \(\tau\). This will allow the precision matrix to have off-diagonal elements different from zero, but will still ensure that it will be positive definite.

# Uncorrelated errors
prior_for_tau         <- prior_diagonal
modelVAR_string_gamma <- glue(modelVAR_string, .open="{{", .close="}}")
jags_output_VAR_gamma <- fit_JAGS(
  modelVAR_string_gamma,
  data_list,
  to_save,
  n_iter=6000,
  n_adapt=1000,
  n_chain=1,
  n_burnin=2000
)

# Possibly correlated errors (Wishart prior)
prior_for_tau           <- prior_wishart
modelVAR_string_wishart <- glue(modelVAR_string, .open="{{", .close="}}")
jags_output_VAR_wishart <- fit_JAGS(
  modelVAR_string_wishart,
  data_list,
  to_save,
  n_iter=6000,
  n_adapt=1000,
  n_chain=1,
  n_burnin=2000
)

# Logging results of the two models
print(glue("Posterior means using diagonal prior"))

## Posterior means using diagonal prior

print(glue("\r\nA"))

## 
## A

print(jags_output_VAR_gamma$mean$A[,,1])

##            [,1]        [,2]
## [1,] 0.91016054 -0.03595962
## [2,] 0.07820191  0.88709689

print(glue("\r\nsigma2"))

## 
## sigma2

print(jags_output_VAR_gamma$mean$sigma2)

##        [,1]      [,2]
## [1,] 2.5992 0.0000000
## [2,] 0.0000 0.9579668

print(glue("\r\n\r\nPosterior means using Wishart prior"))

## 
## 
## Posterior means using Wishart prior

print(glue("\r\nA"))

## 
## A

print(jags_output_VAR_wishart$mean$A[,,1])

##            [,1]       [,2]
## [1,] 0.90728749 -0.0365918
## [2,] 0.07803866  0.8864918

print(glue("\r\nsigma2"))

## 
## sigma2

print(jags_output_VAR_wishart$mean$sigma2)

##           [,1]      [,2]
## [1,] 2.5813639 0.5508038
## [2,] 0.5508038 0.9538458

From these results we see that the noises of the 2 time series are somewhat correlated. We can compute the Pearson Correlation Coefficient to further inspect this hypotheses:

\[ \rho(err_{gdp},err_{cpi}) \;=\;\frac{\text{Cov}(err_{gdp},err_{cpi})}{\sqrt{\text{Var}(err_{gdp})\,\times\,\text{Var}(err_{cpi})}} \;=\; \frac{0.5529124}{\sqrt{2.5870927 \times 0.9520914}} \;\approx\;0.3523 \]

This indicates a moderate positive linear relationship between the two variables. It suggests that as one variable increases, the other tends to increase as well, but not perfectly.

On the other hand, the effects of the 2 different priors on the matrix A is very moderate.

This suggest that the two variables are moderately related through their contemporaneous shocks, but this relation is not dominant enough to significantly affect the lagged relationship between the two variables.

Metrics

At this point, we display all so far computer metrics for all models and both time-series, so as to be able to perform model selection objectively.

# GDP table
kable(
  metrics_df_GDP,
  caption = "Univariate model metrics — GDP series",
  digits  = c(1,1,1,1,3,3,1,1,3),
  format  = "html",
  table.attr = 'cellpadding="10" cellspacing="0"'
) %>%
  kable_styling(full_width = F, position = "center") %>%
  row_spec(0, background = "#000000", color = "#FFFFFF") %>%  # intestazione
  row_spec(1:nrow(metrics_df_GDP), background = "#000000", color = "#FFFFFF")  # righe dati

Univariate model metrics — GDP series

Model	DIC	WAIC	BIC	MSE	MAE	MAPE	SMAPE	MASE
AR(1)	745.0	705.8	691.1	7.985	1.387	75.3	26.8	1.203
AR(2)	686.2	666.8	666.6	12.166	1.826	143.3	31.8	1.583
AR(3)	690.3	666.7	668.9	11.223	1.722	134.3	30.4	1.494
AR(4)	683.4	664.7	672.7	11.807	1.753	133.8	30.0	1.520
AR(5)	646.1	631.7	647.3	8.566	1.534	123.9	26.6	1.330
AR(6)	645.7	631.9	653.4	8.276	1.524	121.9	26.8	1.321
AR(7)	650.7	633.4	657.6	8.240	1.536	121.6	27.0	1.332
AR(8)	656.6	637.2	662.8	8.423	1.581	126.1	28.0	1.371
AR(9)	638.3	622.8	654.7	7.611	1.481	121.6	25.9	1.284
MA(1)	846.7	835.2	838.6	9.168	2.083	114.2	42.7	1.807
MA(2)	791.9	780.3	786.9	15.454	2.483	129.7	44.8	2.153
MA(3)	658.4	638.3	643.8	5.057	1.312	99.5	27.7	1.138
MA(4)	619.8	604.1	612.9	5.911	1.395	134.3	28.1	1.210
MA(5)	653.1	620.4	629.3	7.061	1.393	132.6	25.6	1.208
ARMA(1, 1)	697.7	673.6	672.9	11.412	1.884	136.1	34.8	1.634
ARMA(2, 1)	684.8	664.2	669.9	11.632	1.786	138.6	31.4	1.549
ARMA(3, 1)	676.1	651.2	656.0	11.100	1.781	168.2	32.5	1.544
ARMA(4, 1)	682.7	658.9	665.9	9.473	1.696	140.0	32.3	1.471
ARMA(1, 2)	684.8	675.1	660.2	6.668	1.424	113.0	28.1	1.235
ARMA(2, 2)	646.0	631.2	639.4	7.277	1.550	135.7	30.3	1.344
ARMA(3, 2)	632.8	613.9	626.6	7.647	1.571	132.3	30.3	1.362
ARMA(4, 2)	646.9	624.1	636.8	7.711	1.564	129.2	30.4	1.357
ARMA(1, 3)	591.4	584.6	598.0	6.238	1.291	142.7	24.8	1.119
ARMA(2, 3)	588.2	575.0	588.0	5.544	1.127	133.2	22.2	0.977
ARMA(3, 3)	598.0	579.2	595.2	5.819	1.200	137.5	23.5	1.041
ARMA(4, 3)	612.9	590.1	609.4	5.925	1.249	136.0	24.2	1.084
ARMA(1, 4)	602.7	582.6	592.9	6.152	1.267	139.6	24.4	1.099
ARMA(2, 4)	615.5	584.6	597.1	6.331	1.300	143.6	25.0	1.127
ARMA(3, 4)	608.5	586.0	602.2	6.101	1.255	138.1	24.2	1.088
ARMA(4, 4)	617.4	586.7	606.1	6.198	1.270	138.8	24.3	1.101
VAR_univar_GDP(1)	741.4	704.2	726.5	8.032	1.386	76.9	27.0	1.202

# CPI table
kable(
  metrics_df_CPI,
  caption = "Univariate model metrics — CPI series",
  digits  = c(1,1,1,1,3,3,1,1,3),
  format  = "html",
  table.attr = 'cellpadding="10" cellspacing="0"'
) %>%
  kable_styling(full_width = F, position = "center") %>%
  row_spec(0, background = "#000000", color = "#FFFFFF") %>%
  row_spec(1:nrow(metrics_df_CPI), background = "#000000", color = "#FFFFFF")

Univariate model metrics — CPI series

Model	DIC	WAIC	BIC	MSE	MAE	MAPE	SMAPE	MASE
AR(1)	568.9	552.3	550.4	0.680	0.600	200.0	33.2	0.889
AR(2)	562.2	548.6	551.5	0.579	0.591	212.2	35.8	0.876
AR(3)	560.7	547.5	555.0	0.557	0.576	205.4	37.7	0.854
AR(4)	540.5	532.5	546.6	0.562	0.568	172.4	36.7	0.843
AR(5)	507.1	502.2	519.7	0.500	0.569	149.6	36.1	0.843
AR(6)	507.2	497.9	518.0	0.505	0.575	179.9	36.1	0.852
AR(7)	508.7	497.4	520.5	0.509	0.575	186.2	35.8	0.852
AR(8)	511.9	500.4	528.9	0.514	0.577	198.8	35.9	0.856
AR(9)	481.6	467.0	498.9	0.454	0.532	270.6	34.6	0.789
MA(1)	800.5	792.6	797.1	1.872	1.179	923.2	52.6	1.748
MA(2)	737.7	725.9	726.5	1.316	0.972	646.3	47.2	1.441
MA(3)	571.5	560.7	571.2	0.587	0.591	440.4	33.7	0.876
MA(4)	519.7	508.8	519.4	0.559	0.580	436.5	35.8	0.861
MA(5)	551.9	526.4	536.0	0.474	0.532	331.6	34.9	0.789
ARMA(1, 1)	557.7	546.8	551.1	0.608	0.605	207.9	35.3	0.896
ARMA(2, 1)	549.2	527.8	530.2	0.615	0.627	266.7	41.0	0.929
ARMA(3, 1)	545.7	540.9	556.5	0.557	0.571	194.1	37.0	0.846
ARMA(4, 1)	517.2	500.5	509.7	0.542	0.614	270.3	40.8	0.910
ARMA(1, 2)	543.2	533.6	533.2	0.535	0.531	151.4	34.9	0.787
ARMA(2, 2)	547.3	526.6	534.4	0.645	0.644	302.4	41.8	0.955
ARMA(3, 2)	501.4	494.5	508.9	0.471	0.508	139.9	35.2	0.754
ARMA(4, 2)	517.6	516.2	514.0	0.433	0.517	176.7	33.9	0.767
ARMA(1, 3)	475.8	467.4	477.5	0.424	0.535	301.9	34.9	0.793
ARMA(2, 3)	470.0	457.2	471.7	0.400	0.520	259.7	34.9	0.772
ARMA(3, 3)	433.3	421.3	438.7	0.379	0.495	236.4	34.6	0.735
ARMA(4, 3)	468.4	448.5	466.9	0.434	0.526	265.8	36.2	0.779
ARMA(1, 4)	466.7	450.1	460.4	0.400	0.491	229.6	33.4	0.728
ARMA(2, 4)	475.2	448.2	457.7	0.437	0.495	222.0	34.0	0.735
ARMA(3, 4)	467.4	453.5	471.1	0.434	0.493	220.5	33.8	0.731
ARMA(4, 4)	460.0	436.4	457.6	0.430	0.506	203.3	34.9	0.750
VAR_univar_GDP(1)	548.9	537.2	567.6	0.720	0.638	204.1	39.6	0.945

NOTICE: The MAPE metric exceeds 100% because the series (quarterly or monthly GDP growth rates) has true values (denominators) that are very small or close to zero. This results in an amplified percentage.

Trace plots and posterior densities

We will now display some trace plots for some parameters of some of the fitted models and relative posterior density plots to make sure the MCMC samples are not getting stuck in some local regions, and are behaving well overall.

AR(9) trace and density plots (GDP)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_GDP,
  "order"=9,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelAR_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 251
##    Unobserved stochastic nodes: 352
##    Total graph size: 1901
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:11])

mcmc_dens(output$samples[,1:11])

mcmc_intervals(output$samples[,1:11])

MA(3) trace and density plots (GDP)

n_chain = if (intense) 3 else 1

to_save <- c(
  "beta",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_GDP,
  "order"=3,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelMA_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 257
##    Unobserved stochastic nodes: 397
##    Total graph size: 2265
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

mcmc_trace(output$samples[,1:5])

mcmc_dens(output$samples[,1:5])

mcmc_intervals(output$samples[,1:5])

ARMA order selection and traces (GDP)

## --- 1. Support Function: matrix [p × q] ---------------------------
make_matrix <- function(df, metric, P, Q){
  M <- matrix(NA_real_, nrow = P, ncol = Q)   # row = p,  column = q
  for(i in seq_len(nrow(df))){
    nums <- as.integer(regmatches(df$Model[i],
                   gregexpr("\\d+", df$Model[i]))[[1]])
    if(length(nums) == 2){
      p <- nums[1];  q <- nums[2]
      M[p, q] <- df[i, metric]
    }
  }
  M
}
## --- 2. function for plotting axis  -------------------------
plot_heat <- function(mat, title){
  image(x = 1:ncol(mat),            
        y = 1:nrow(mat),            
        z = mat,                  
        col = colorRampPalette(
                c("navy","deepskyblue","yellow","orange","red"))(100),
        axes = FALSE, xlab = "", ylab = "", main = title, useRaster = TRUE)
  contour(1:ncol(mat), 1:nrow(mat), mat,
          add = TRUE, lwd = 0.7, drawlabels = FALSE)
  axis(1, at = 1:ncol(mat), labels = 1:ncol(mat), cex.axis = 0.75) # q
  axis(2, at = 1:nrow(mat), labels = 1:nrow(mat), las = 1,
       cex.axis = 0.75)                                           # p
  box()
}

## --- 3. parameters & palette --------------------------------------------
max_p <- max_order_ar   # max order AR
max_q <- max_order_ma   # max order MA
metrics <- c("DIC","WAIC","BIC","MSE",
             "MAE","MAPE","SMAPE","MASE")

## --- 4. heat-map 2×4 for axis p↔q --------------------------------------
op <- par(mfrow = c(2,4), mar = c(2.2,2.2,2.4,0.6))

for(met in metrics){
  mat <- make_matrix(metrics_df_GDP, met, max_p, max_q)
  plot_heat(mat, paste("ARMA –", met))
}

par(op)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2"
)

order_ar <- 2
order_ma <- 3
data_list <- list(
  "y"=series_GDP,
  "N_train"=N_train,
  "N_test"=N_test,
  "order"=max(order_ar, order_ma),
  "order_ar"=order_ar,
  "order_ma"=order_ma,
  "train_metrics_skip"=50
)

# Perform simulations and visualise output
output <- fit_JAGS(
  modelARMA_string,
  data_list,
  to_save,
  n_thin=50,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=50000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 257
##    Unobserved stochastic nodes: 399
##    Total graph size: 2964
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 49000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:8])

mcmc_dens(output$samples[,1:8])

mcmc_intervals(output$samples[,1:8])

AR(9) trace and density plots (CPI)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_CPI,
  "order"=9,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelAR_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 251
##    Unobserved stochastic nodes: 352
##    Total graph size: 1901
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:11])

mcmc_dens(output$samples[,1:11])

mcmc_intervals(output$samples[,1:11])

MA(4) trace and density plots (CPI)

n_chain = if (intense) 3 else 1

to_save <- c(
  "beta",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_CPI,
  "order"=4,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelMA_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 256
##    Unobserved stochastic nodes: 397
##    Total graph size: 2261
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

mcmc_trace(output$samples[,1:6])

mcmc_dens(output$samples[,1:6])

mcmc_intervals(output$samples[,1:6])

ARMA order selection and traces (CPI)

## --- 1. Support Function: matrix [p × q] ---------------------------
make_matrix <- function(df, metric, P, Q){
  M <- matrix(NA_real_, nrow = P, ncol = Q)  
  for(i in seq_len(nrow(df))){
    nums <- as.integer(regmatches(df$Model[i],
                   gregexpr("\\d+", df$Model[i]))[[1]])
    if(length(nums) == 2){
      p <- nums[1];  q <- nums[2]
      M[p, q] <- df[i, metric]
    }
  }
  M
}

## --- 2. function for plotting axis  -------------------------
plot_heat <- function(mat, title){
  image(x = 1:ncol(mat),         
        y = 1:nrow(mat),          
        z = mat,                   
        col = colorRampPalette(
                c("navy","deepskyblue","yellow","orange","red"))(100),
        axes = FALSE, xlab = "", ylab = "", main = title, useRaster = TRUE)
  contour(1:ncol(mat), 1:nrow(mat), mat,
          add = TRUE, lwd = 0.7, drawlabels = FALSE)
  axis(1, at = 1:ncol(mat), labels = 1:ncol(mat), cex.axis = 0.75) # q
  axis(2, at = 1:nrow(mat), labels = 1:nrow(mat), las = 1,
       cex.axis = 0.75)                                           # p
  box()
}

## --- 3. parameters & palette --------------------------------------------
max_p <- max_order_ar   # max order AR
max_q <- max_order_ma   # max order MA
metrics <- c("DIC","WAIC","BIC","MSE",
             "MAE","MAPE","SMAPE","MASE")

## --- 4. heat-map 2×4 for axis p↔q --------------------------------------
op <- par(mfrow = c(2,4), mar = c(2.2,2.2,2.4,0.6))

for(met in metrics){
  mat <- make_matrix(metrics_df_CPI, met, max_p, max_q)
  plot_heat(mat, paste("ARMA –", met))
}

par(op)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2"
)

order_ar <- 3
order_ma <- 3
data_list <- list(
  "y"=series_CPI,
  "N_train"=N_train,
  "N_test"=N_test,
  "order"=max(order_ar, order_ma),
  "order_ar"=order_ar,
  "order_ma"=order_ma,
  "train_metrics_skip"=50
)

# Perform simulations and visualise output
output <- fit_JAGS(
  modelARMA_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 257
##    Unobserved stochastic nodes: 400
##    Total graph size: 2965
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:8])

mcmc_dens(output$samples[,1:8])

mcmc_intervals(output$samples[,1:8])

VAR(1) trace and density plots

to_save <- c(
  "A",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=y_var,
  "order"=1,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  final_modelVAR_string,
  data_list,
  to_save,
  n_thin=3,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=6000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 259
##    Unobserved stochastic nodes: 356
##    Total graph size: 2735
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 5000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

mcmc_trace(output$samples[,1:10])

mcmc_dens(output$samples[,1:10])

mcmc_intervals(output$samples[,1:10])

Model selection

Based on the calculated metrics, the models that seem to perform best are \(\text{ARMA}(2, 3)\) and \(\text{ARMA}(3, 3)\), both for the GDP and inflation time series.

An honorable mention must be given to the \(\text{AR}(9)\) model for the CPI time series, where the out-of-sample prediction for the next year and a half, given only the in-sample data, is almost perfectly aligned with the actual values, as shown in the following plot.

if (intense) {
  visualise(x_axis, series_CPI, output_ar9_cpi, title=glue("CPI series, AR({order})"))
}

Anomaly detection

In this last section we will examine anomalies and outliers which are present in our two time series. So far, we performed model selection and predictions of future values for the series, however we could clearly see that both the training and test set presented samples which deviated quite heavily from their neighboring values. This can be cause of disturbances during the learning phase and also a reduced accuracy during forecasting.

The objective of this section is to detect said outliers and visualise them so as to clearly show which samples may be problematic during fitting and out-of-sample prediction.

In order to do this, for each timestamp \(t\), we consider the incremental values of our series, each of which is to be treated as an independent random variable:

\[ \Delta_t := y_{t} - y_{t-1} \stackrel{iid}{\sim} \mathcal{N}(\mu, \sigma_t^2) \]

In order to be able to perform anomaly detection, we must ensure that the sequence of incremental values \(\Delta_t\) can be approximated to a random walk. This can be verified by trying to fit an \(\text{AR}(1)\) process to the initial sequence of values \(y_t\) and checking whether the value \(|\alpha|\) is close to 1:

# Data lists for both delta-series
data_list_GDP <- list(
  "y"=series_GDP,
  "train_metrics_skip"=50,
  "order"=1,
  "N_train"=N_tot,
  "N_test"=1          # No need to perform forecasting
)
data_list_CPI <- list(
  "y"=series_CPI,
  "train_metrics_skip"=50,
  "order"=1,
  "N_train"=N_tot,
  "N_test"=1          # No need to perform forecasting
)

# Only need to inspect value of alpha for both series
to_save <- c(
  "alpha"
)

# Fit AR(1) models
output_GDP <- fit_JAGS(
  modelAR_string,
  data_list_GDP,
  to_save
)
output_CPI <- fit_JAGS(
  modelAR_string,
  data_list_CPI,
  to_save
)

# Inspecting values of alpha
print(glue("Alpha value for GDP has posterior mean {output_GDP$mean$alpha} with standard deviation of {output_GDP$sd$alpha}"))

## Alpha value for GDP has posterior mean 0.85760849050055 with standard deviation of 0.0304596430365325

print(glue("Alpha value for CPI has posterior mean {output_CPI$mean$alpha} with standard deviation of {output_CPI$sd$alpha}"))

## Alpha value for CPI has posterior mean 0.938589495371746 with standard deviation of 0.0197994927970986

Now that we know that the two series can be approximated reasonably well to random walks (\(\alpha_{GDP} = 0.86\) and \(\alpha_{CPI} = 0.938\)) we proceede with our anomaly detection.

Notice that the random variable \(\Delta_t\) has a fixed mean value \(\mu\), but a time dependant variance \(\sigma_t^2\), which we model as:

\[ \sigma_t^{-2} = \tau_t = \beta_1 + \beta_2 \delta_t \qquad \delta_t \sim Ber(p) \]

We can, therefore, define the following JAGS model-string which will fit the individual values of \(\delta_t\).

anomaly_detection_string <- "
model {

  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################

  for(t in 1:N_tot) {
    delta[t]  ~  dnorm(mu, tau[t])
    tau[t]    <- beta[1] + beta[2] * d[t]
  }
  
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  # Beta coefficients, Bernoulli probability and mean
  beta[1] ~ dexp(0.1)
  beta[2] ~ dexp(0.1)
  p       ~ dbeta(1, 1)
  mu      ~ dnorm(0, 0.01)
  
  # Delta value
  for(t in 1:N_tot){
    d[t] ~ dbern(p)
  }
}"

# Obtaining series of differences
delta_GDP <- series_GDP[2:N_tot] - series_GDP[1:(N_tot-1)]
delta_CPI <- series_CPI[2:N_tot] - series_CPI[1:(N_tot-1)]

# Setting up parameters for JAGS
data_list_GDP <- list(
  "delta"=delta_GDP,
  "N_tot"=N_tot-1
)
data_list_CPI <- list(
  "delta"=delta_CPI,
  "N_tot"=N_tot-1
)
to_save <- c(
  "d",
  "mu",
  "tau"
)

# Perform simulations
anomaly_output_GDP <- fit_JAGS(
  anomaly_detection_string,
  data_list_GDP,
  to_save,
  n_iter=10000
)
anomaly_output_CPI <- fit_JAGS(
  anomaly_detection_string,
  data_list_CPI,
  to_save,
  n_iter=10000
)

Once obtained both \(\{\delta_t\}_{t=1...N}\) sequences we can classify their values according to their expected values. Concretely, we classify a timestamp \(t\) such that \(\mathbb{E}[\delta_t] < 0.2\) as “high variance” (where strong oscillations may occur) and timestamps \(t\) with \(0.2 < \mathbb{E}[\delta_t] < 0.6\) as samples with “medium variance” (where we still might have an anomaly, but of weaker nature). All remaining samples can be considered stable and only subject to natural stochastic oscillations.

# Classifying anomalies for both series
anomalies_GDP <- classify_anomalies(anomaly_output_GDP, N_tot, title="GDP")

## Found 36 strong and 27 medium outliers out 305 total samples for GDP

anomalies_CPI <- classify_anomalies(anomaly_output_CPI, N_tot, title="CPI")

## Found 23 strong and 42 medium outliers out 305 total samples for CPI

# Displaying graphically the classification
delta_GDP <- anomaly_output_GDP$mean$d
delta_CPI <- anomaly_output_CPI$mean$d

display_classification(
  x_axis,
  delta_GDP,
  anomalies_GDP$high_var_indices,
  anomalies_GDP$med_var_indices,
  anomalies_GDP$low_var_indices,
  title="GDP"
)

display_classification(
  x_axis,
  delta_CPI,
  anomalies_CPI$high_var_indices,
  anomalies_CPI$med_var_indices,
  anomalies_CPI$low_var_indices,
  title="CPI"
)

Finally, we can plot both our time series displaying the points in time where high variances were observed: red vertical lines display timestamps where strong anomalies were detected, while orange dotted lines show medium anomalies.

# Zoom value: we only plot starting from this point on so as not to clog up the graph
start <- 150

visualise_anomalies(
  x_axis,
  series_GDP,
  anomalies_GDP$high_var_indices,
  anomalies_GDP$med_var_indices,
  title="GDP",
  start=start
)

visualise_anomalies(
  x_axis,
  series_CPI,
  anomalies_CPI$high_var_indices,
  anomalies_CPI$med_var_indices,
  title="CPI",
  start=start
)